/*
    KMC Simulation
    Copyright (c) 2012 Filip Rozbořil

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "includes.h"
#include "constants.h"
#include "KmcSimulation.h"
#include "MersenneTwister.h"
#include "Grid.h"
#include "Utils.h"
#include "PartialEvent.h"
#include "HopEvent.h"
#include "AtomDepositionEvent.h"
#include "Histogram.h"

const uint32_t KmcSimulation :: Version = UINT32_C(4000);

KmcSimulation :: KmcSimulation(KmcSettings &sett) : mGrid(sett), mTime(0.), mTotalEvents(UINT64_C(0)), mRateUpdates(UINT64_C(0)), mDepositionHasEnded(false) 
{
	uint32_t seed;
	sett.getValue("RandomSeed", seed);
	mRand = new MersenneTwister(seed);
	
	sett.getValue("DepositionTime", mDepositionTime);
	assert(mDepositionTime > 0.0);
	
	sett.getValue("RelaxationTime", mRelaxationTime);
	assert(mRelaxationTime >= 0.0);
	
	sett.getValue("MaxEvents", mMaxEvents);
	assert(mMaxEvents > 0);
	
	mTotalTime = mDepositionTime + mRelaxationTime;
	
	mProxyEvents.reserve(Site::Count);
	mEvents.reserve(Site::Count + 1);

	for(Site::type from = 0; from < Site::Count; from++)
 	{
 		auto fromEvent = new PartialEvent(Site::Count);
 		for(Site::type to = 0; to < Site::Count; to++)
 		{
 			auto toEvent = new PartialEvent(Site::Count);
 			toEvent->addSubEvent(new HopEvent<Direction::Down>(from, to, mGrid));
 			toEvent->addSubEvent(new HopEvent<Direction::Up>(from, to, mGrid));
 			toEvent->addSubEvent(new HopEvent<Direction::Left>(from, to, mGrid));
 			toEvent->addSubEvent(new HopEvent<Direction::Right>(from, to, mGrid));
 			fromEvent->addSubEvent(toEvent);
 		}
 		mProxyEvents.push_back(fromEvent);
 		mEvents.push_back(fromEvent);	
 	}

 	mDeposition = new AtomDepositionEvent(mGrid.getCoverage() / mDepositionTime * mGrid.getElements(), mGrid);
 	mEvents.push_back(mDeposition);
 	
 	mTotalRate = mDeposition->getRate();

}

KmcSimulation :: ~KmcSimulation()
{ 
	delete mRand;
	for(auto e : mEvents)
	{
		delete e;
	}
}

void KmcSimulation :: updateRates()
{
	Trace("Updating rates");
	mRateUpdates++;
	mTotalRate = 0.;
	for(auto e : mProxyEvents)
	{
		e->updateRate();
		mTotalRate += e->getRate();
	}
	mTotalRate += mDeposition->getRate();
}

void KmcSimulation :: run()
{
	cout << "[Runtime details]" << endl;
	cout << mSw.getElapsed(0.) << " Simulation started." << endl;
	
	const int totalSteps = 10;
	uint64_t maxEvents = mMaxEvents;
	uint64_t eventStep = mMaxEvents / totalSteps;
	double eventTimeStep = mTotalTime / totalSteps;
	double timeStep = mDepositionTime / totalSteps;
	
	//start stopwatch
	mSw.start();
	
	// The evolution of the system is not important before there are at least 2 atoms on the grid.
	// Between the first and the second atom deposition, there are a lot of unnecessary hop events that
	// just randomize the position of the first atom on an empty grid.
	// We improve performance by starting with 2 atoms at random positions.
	
	//deposit two adatoms
	mDeposition->handle(mRand);
	mTime += -log(mRand->getDouble())/mTotalRate;
	mDeposition->handle(mRand);
	mTime += -log(mRand->getDouble())/mTotalRate;

	Trace("Initial time after two atom depositions: " << mTime);

 	while(mTime < mTotalTime && mTotalEvents < mMaxEvents)
 	{				
 	 	if (!mDepositionHasEnded && (mTime > mDepositionTime))
 	 	{        
			//make the number of atoms an even number, so they can all dimerize without prolonging the simulation
			if((mGrid.getAtomCount() % 2) != 0)
			{ 
				mDeposition->handle(mRand);
			}
			mDeposition->setRate(0.);
			mDepositionHasEnded = true;
			updateRates();
			
			cout << "Deposition Finished. Time: " << mTime << "/" << mTotalTime << endl;
			
			//separate event percentage count for relaxation
			maxEvents -= mTotalEvents;
			eventStep = mTotalEvents + (totalSteps - mTotalEvents % totalSteps) + maxEvents / totalSteps;
		}
		
		if(mGrid.hasStateChanged())
		{
			Trace("[KMC] Grid state has changed to: " << mGrid.getHash());
			
			updateRates();
			
			//Trace("[KMC] Event rates:");
			Trace("[KMC] Total rate: " << mTotalRate);
		}		

		double deltaTime = -log(mRand->getDouble())/mTotalRate;
		mTime += deltaTime;
		
		Trace("[KMC] Random delta time: " << deltaTime);
		
		if(mTime > mTotalTime)
		{
		 	break;
		}
		
 		chooseEvent()->handle(mRand);
 		
 		Trace("[KMC] Total events: " << mTotalEvents);
 		
#ifdef GRIDCHECK
 		mGrid.checkConsistency();
#endif
 		
 		if(mTime >= timeStep)
 		{
 			if(mDepositionHasEnded)
 			{
 				double percent = (timeStep - mDepositionTime) / mRelaxationTime;
 				mSw.mark();
 				double seconds = mSw.getSeconds();
 				cout << mSw.getElapsed(seconds) << " " << 100 * percent << "% of relaxation (" << mTotalEvents / seconds << " ev/s)" << endl;
 				timeStep += mRelaxationTime / totalSteps;
 				eventStep += maxEvents / totalSteps;
 			}
 			else
 			{
 			 	double percent = timeStep / mDepositionTime;
 			 	mSw.mark();
 				double seconds = mSw.getSeconds();
 				cout << mSw.getElapsed(seconds) << " " << 100 * percent << "% of deposition (" << mTotalEvents / seconds << " ev/s)" << endl;
 				timeStep += (percent >= 1.0 ? mRelaxationTime : mDepositionTime) / totalSteps;
 				eventStep += maxEvents / totalSteps;
 			}
 		}
 		if(mTime >= eventTimeStep)
 		{
 			eventTimeStep += mTotalTime / totalSteps;
 			eventStep += maxEvents / totalSteps;
 		}
 		if(mTotalEvents >= eventStep)
 		{
 			double percent = (double)eventStep / mMaxEvents;
 			mSw.mark();
 			double seconds = mSw.getSeconds();
 			cout << mSw.getElapsed(seconds) << " " << 100 * percent << "% of events (" << mTotalEvents / seconds << " ev/s), Time: " << mTime << endl;
 			timeStep += (mDepositionHasEnded ? mRelaxationTime : mDepositionTime) / totalSteps;
 			eventStep += maxEvents / totalSteps;
 			eventTimeStep += mTotalTime / totalSteps;
 		}
 		
 	}
 	mSw.mark();
	cout << " Simulation finished."; 
	outputDetails(cout);
}

inline Event* KmcSimulation :: chooseEvent()
{
	double pivot = mRand->getDouble() * mTotalRate;
	Trace("[KMC] Choosing event with cumulative rate " << pivot);
	
	auto i = mEvents.begin();
 	double probe = (*i)->getRate();
 	
	while(probe < pivot)
	{
		++i;
		probe += (*i)->getRate();
	}
	
	mTotalEvents++;
	return (*i);
}

void KmcSimulation :: outputDetails(ostream& os)
{
	double seconds = mSw.getSeconds();
 	os 	<< "Elapsed time: " << mSw.getElapsed(seconds) << endl
 		<< "Simulation time: " << mTime << "/" << mTotalTime << endl
 		<< "Total events: " << mTotalEvents << endl
 		<< "Events per second: " << mTotalEvents / seconds << endl;
 		
 	for(auto e : mEvents)
	{
	 	e->printCount(cout);
	}
	
	Histogram hist(mGrid);
	
	os 	<< "Total rate updates: " << mRateUpdates << endl 
	 	<< endl
	 	<< "[Simulation results]" << endl 
		<< "Real coverage: " << setprecision(3) << ((double)mGrid.getAtomCount())/mGrid.getElements() << " ML" << endl
		<< "Atoms: " << mGrid.getAtomCount() << endl
		<< "Chains: " << hist.getChainCount() << endl
		<< "Kinks: " << mGrid.getKinkCount() << endl
		<< endl
		<< "[Histogram]" << endl
		<< hist << endl
		<< "[Grid]" << endl
		<< mGrid; 

	os << flush;
}
